
clear
set more off
program drop _all

local owndir = "" // put your path here

global dirdata = "`owndir'/Raw_Data"
global dirin = "`owndir'/Dta"
global dirout = "`owndir'/Output"



************
* PI_ij calculation
************

*** Step 1: Share of people working at home by district before pandemic

use "$dirin/Seoul_Survey_clean.dta",clear
drop if de4==9 // no income

gen age_group=0
replace age_group=1 if de2>=2&de2<=5 // 20-30, 30-40, 40-50, 50-60
replace age_group=2 if de2==6 // 60+
drop if age_group==0 // 10-20

gen work_at_home=0
replace work_at_home=1 if destination=="Not commuting"
collapse (sum) wtb1, by (work_at_home district age_group)
bys district age_group: egen sum_wtb1=sum(wtb1)
gen share=wtb1/sum_wtb1
keep if work_at_home==1
keep district age_group share
sort age_group district

save "$dirin/initial_home_share_by_district.dta",replace


*** Step 2: Calculate decline in flow over time by age group

use "$dirin/flows_ij_by_age.dta", clear

collapse (sum) flow, by (all_dates origin age_group)

gen dow=dow(all_dates)
gen month=month(all_dates)
gen day=day(all_dates)

gen flow2=flow if day==2 & month==1
gen flow3=flow if day==3 & month==1
gen flow4=flow if day==4 & month==1
gen flow5=flow if day==5 & month==1
gen flow6=flow if day==6 & month==1
gen flow7=flow if day==7 & month==1
gen flow8=flow if day==8 & month==1

bysort origin age_group: egen min_dow4=min(flow2)
bysort origin age_group: egen min_dow5=min(flow3)
bysort origin age_group: egen min_dow6=min(flow4)
bysort origin age_group: egen min_dow0=min(flow5)
bysort origin age_group: egen min_dow1=min(flow6)
bysort origin age_group: egen min_dow2=min(flow7)
bysort origin age_group: egen min_dow3=min(flow8)

gen normalization=min_dow3 if dow==3
replace normalization=min_dow4 if dow==4
replace normalization=min_dow5 if dow==5
replace normalization=min_dow6 if dow==6
replace normalization=min_dow0 if dow==0
replace normalization=min_dow1 if dow==1
replace normalization=min_dow2 if dow==2

keep origin age_group flow all_dates dow normalization
gen flow_norm=flow/normalization

save "$dirin/flows_from_origin_normalized_by_age_group.dta",replace


*** Step 3: Calculating pi_ij^a

use "$dirin/flows_ij_by_age.dta",clear
gen month=month(all_dates)
rename origin district
merge m:1 district age_group using "$dirin/initial_home_share_by_district.dta", keepusing(share)
drop _merge
rename district origin
gen one_minus_home=1-share
drop share
merge m:1 origin age_group all_dates using "$dirin/flows_from_origin_normalized_by_age_group.dta", keepusing(flow_norm)
drop _merge
replace one_minus_home=one_minus_home*flow_norm
drop flow_norm

bys origin all_dates age_group: egen denominator=sum(flow) // summing 25 destinations
gen pi_ij=flow/denominator*one_minus_home // sum pi_ij over j (=one_minus_home) + pi_ii (1-one_minus_home) = 1

keep all_dates age_group origin destination pi_ij
sort all_dates age_group origin destination
save "$dirin/pi_ij_ii.dta",replace


************
* Gravity Equation Regression
************

* 1. Clean distance

use "$dirdata/d_GIS.dta",clear
reshape long var, i(district_num) j(j)
rename district_num d
merge m:1 d using "$dirdata/district_list.dta",keepusing(district)
drop _merge d
rename district origin
rename j d
merge m:1 d using "$dirdata/district_list.dta",keepusing(district)
drop _merge d
rename district destination
rename var d
save "$dirin/d_GIS_clean.dta",replace


* 2. bring distance, case (for past 14 days), and visit (for past 14 days)

use "$dirin/pi_ij_ii.dta",clear
merge m:1 origin destination using "$dirin/d_GIS_clean.dta", keepusing(d)
drop _merge
gen month=month(all_dates)
gen day=day(all_dates)
rename destination district

merge m:1 month day using "$dirin/covid19_Seoul_cum_cases_14.dta", keepusing(cum_cases_total_14)
drop if _merge==2
drop _merge
replace cum_cases_total_14=0 if cum_cases_total_14==.
merge m:1 district month day using "$dirin/covid19_Seoul_cum_cases_14_by_district.dta",keepusing(cum_cases_14)
drop if _merge==2
drop _merge
replace cum_cases_14=0 if cum_cases_14==.
merge m:1 district month day using "$dirin/covid19_Seoul_cum_visits_14_by_district.dta",keepusing(cum_visits_14)
drop if _merge==2
drop _merge
replace cum_visits_14=0 if cum_visits_14==.
rename district destination

gen ln_pi=ln(pi_ij)
gen ln_cum_cases_14=ln(cum_cases_14+1)
gen ln_cum_cases_14_squared=ln_cum_cases_14*ln_cum_cases_14
gen ln_cum_visits_14=ln(cum_visits_14+1)
gen ln_cum_visits_14_squared=ln_cum_visits_14*ln_cum_visits_14
gen ln_cum_cases_total_14=ln(cum_cases_total_14+1)
gen ln_cum_cases_total_14_squared=ln_cum_cases_total_14*ln_cum_cases_total_14
gen net_visits_14= cum_visits_14 - cum_cases_14 
replace net_visits_14=0 if net_visits_14<0
gen ln_net_visits_14=ln(net_visits_14 + 1)
gen ln_net_visits_14_squared=ln_net_visits_14*ln_net_visits_14

egen origin_group=group(origin)
egen destination_group=group(destination)
egen date_group=group(all_dates)
gen  dow=dow(all_dates)
gen weekend=(dow==6|dow==0)

* 3. run kappa, delta and ksi, age-specific (not weekday/weekend specific), linear

gen ln_cum_cases_14xWeekend=ln_cum_cases_14*weekend
gen ln_net_visits_14xWeekend=ln_net_visits_14*weekend
gen ln_cum_cases_total_14xWeekend=ln_cum_cases_total_14*weekend

// epsilon_weekday=4.1642;
// epsilon_weekend=4.9144;

constraint 1 ln_cum_cases_14xWeekend = (4.9144-4.1642)/4.1642 * ln_cum_cases_14
constraint 2 ln_net_visits_14xWeekend = (4.9144-4.1642)/4.1642 * ln_net_visits_14
constraint 3 ln_cum_cases_total_14xWeekend = (4.9144-4.1642)/4.1642 * ln_cum_cases_total_14

// epsilon^k delta  X + (epsilon^d - epsilon^k) delta X * weekend
// beta /epsilon^k = delta 

preserve
keep if month==1&day==1
rename pi_ij base_pi_ij
keep origin destination age_group base_pi_ij
save "$dirin/base_pi.dta",replace
restore

merge m:1 origin destination age_group using "$dirin/base_pi.dta"
gen diff_ln_pi=ln(pi/base_pi_ij)

cnsreg diff_ln_pi ln_cum_cases_14 ln_cum_cases_14xWeekend ln_net_visits_14 ln_net_visits_14xWeekend weekend  i.destination_group i.all_dates if age_group==1, constraint(1,2) vce(bootstrap, reps(50) seed(1000) cluster(all_dates destination_group))
// delta^young=.0087075/4.1642=0.00209103789, ksi^young=.0057727/4.1642=0.00138626867



cnsreg diff_ln_pi ln_cum_cases_14 ln_cum_cases_14xWeekend ln_net_visits_14 ln_net_visits_14xWeekend weekend  i.destination_group i.all_dates if age_group==2, constraint(1,2) vce(bootstrap, reps(50) seed(1000) cluster(all_dates destination_group))
// delta^old=.0103115/4.1642=0.00247622592, ksi^old=.0040247/4.1642=0.00096650016


* saving time FEs
preserve
cnsreg diff_ln_pi ln_cum_cases_14 ln_cum_cases_14xWeekend ln_net_visits_14 ln_net_visits_14xWeekend weekend  i.destination_group i.all_dates if age_group==1, constraint(1,2)
	matrix results=[e(b)',vecdiag(e(V))']
	svmat results
	rename results1 coef
	rename results2 se
	keep coef se
	drop if coef == .
	gen dum=_n
	drop if dum<32
	drop if dum>183
	drop dum
	gen coef_u = coef + 1.96*sqrt(se)
	gen coef_l = coef - 1.96*sqrt(se)
	gen d=_n
	save "$dirin/timeFE_young.dta",replace
restore

preserve
cnsreg diff_ln_pi ln_cum_cases_14 ln_cum_cases_14xWeekend ln_net_visits_14 ln_net_visits_14xWeekend weekend  i.destination_group i.all_dates if age_group==2, constraint(1,2)
	matrix results=[e(b)',vecdiag(e(V))']
	svmat results
	rename results1 coef
	rename results2 se
	keep coef se
	drop if coef == .
	gen dum=_n
	drop if dum<32
	drop if dum>183
	drop dum
	gen coef_u = coef + 1.96*sqrt(se)
	gen coef_l = coef - 1.96*sqrt(se)
	gen d=_n
	save "$dirin/timeFE_old.dta",replace
restore

use "$dirin/covid19_Seoul_cum_cases_14.dta",clear
gen d=_n+23
merge 1:1 d using "$dirin/timeFE_young.dta", keepusing(coef)
drop _merge
rename coef coef_young
merge 1:1 d using "$dirin/timeFE_old.dta", keepusing(coef)
drop _merge
rename coef coef_old
gen ln_cum_cases_total_14=ln(cum_cases_total_14+1)
reg coef_young ln_cum_cases_total_14 // .0243171/4.1642=0.00583956102
reg coef_old ln_cum_cases_total_14   // .0306678/4.1642=0.00736463186
